load("agg_day_year_all.rdata")
load("n_day_year_all.rdata")

## commented out for replication files 
#  install.packages("xtable")
#  install.packages("dplyr")
#  install.packages("rmeta")
#  install.packages("rdrobust")

library(xtable)
library(dplyr)
library(rmeta)
library(rdrobust)

year <- c(2009, 2013, 2014, 2015)

data <- left_join(data_day_year, data_n_year,
                  by = c("days", "year"))

data_close <- 
  data %>%
  mutate(voted     = round(par_vote*n),
         abstained = round((1-par_vote)*n),
         treated   = days > 0) %>% 
  select(year, days, treated, voted, abstained)

data_voted <- 
  as.data.frame(lapply(data_close, 
                       function(x,p) rep(x,p), 
                       data_close[["voted"]])) %>%
  mutate(voted = 1)

data_abstained <- 
  as.data.frame(lapply(data_close, 
                       function(x,p) rep(x,p), 
                       data_close[["abstained"]])) %>%
  mutate(voted = 0)

data <- 
  rbind(data_voted, data_abstained)

# array with coefficients for years and models 

main_table_out <- matrix(NA, nrow = 11, ncol = 5)

coef_matrix <- matrix(NA, nrow = 4, ncol = 5)
for (i in 1:4){
  year_vec <- year[i]
  data_year <- 
    data %>%
    ungroup() %>%
    mutate(days = days) %>% 
    filter(year == year_vec )
  
  rd <- rdrobust(data_year$voted, data_year$days)
  coef_matrix[c(1,3), i] <- rd$coef[c(1,3)]
  coef_matrix[c(2,4), i] <- rd$se[c(1,3)]
  
  main_table_out[7, i] <- as.numeric(data_year %>% 
                                       filter(treated == 0 & abs(days) < rd$h_l) %>% 
                                       summarise(mean(voted)))
  main_table_out[8  ,i] <- rd$N_h_l
  main_table_out[9  ,i] <- rd$N_h_r
  main_table_out[10 ,i] <- round(rd$h_l)
  main_table_out[11 ,i] <- round(rd$b_l)
  print(i)
}

# Pool coefficients using fixed effects meta analysis
coef_matrix[1:2, 5] <- c(meta.summaries(coef_matrix[ 1, 1:4],
                                        coef_matrix[ 2, 1:4])$summary,
                         meta.summaries(coef_matrix[ 1, 1:4],
                                        coef_matrix[ 2, 1:4])$se.summary)
coef_matrix[3:4, 5] <- c(meta.summaries(coef_matrix[ 3, 1:4],
                                        coef_matrix[ 4, 1:4])$summary,
                         meta.summaries(coef_matrix[ 3, 1:4],
                                        coef_matrix[ 4, 1:4])$se.summary)

main_table_out[7,5] <- weighted.mean(x = main_table_out[7, 1:4],
                                    w = main_table_out[8, 1:4])
main_table_out[7,] <- round(100*main_table_out[7,], 2)
main_table_out[c(8,9),5] <- rowSums(main_table_out[c(8,9), -5])
for (j in 1:5){
  main_table_out[1,j] <- 100*round(coef_matrix[1,j], 4)
  main_table_out[2,j] <- 100*round(coef_matrix[2,j], 4)
  main_table_out[4,j] <- 100*round(coef_matrix[3,j], 4)
  main_table_out[5,j] <- 100*round(coef_matrix[4,j], 4)
  main_table_out[3     , j] <- paste("[",
                                     100*round(coef_matrix[1,j] + coef_matrix[2,j] * qnorm(0.025), 4),
                                     ", ",
                                     100*round(coef_matrix[1,j] + coef_matrix[2,j] * qnorm(0.975), 4),
                                     "]",
                                     sep = "")
  main_table_out[6     , j] <- paste("[",
                                     100*round(coef_matrix[3,j] + coef_matrix[4,j] * qnorm(0.025), 4),
                                     ", ",
                                     100*round(coef_matrix[3,j] + coef_matrix[4,j] * qnorm(0.975), 4),
                                     "]",
                                     sep = "")  
  
}

rownames(main_table_out) <- c("Estimate", "St.err.", "CI",
                              "Robust estimate", "Robust st.err.", "Robust CI",
                              "Mean turnout among parents of ineligible", "N_ineligible",
                              "N_eligible", "Bandwidth", "Bias-correction bandwidth")
colnames(main_table_out) <- c("2009", "2013", "2014", "2015", "Pooled Effect")

main_table_out_all <- main_table_out

# Repeat for only 32 recurrent municipalities

load("agg_day_year_allyears.rdata")
load("n_day_year_allyears.rdata")
#rename dataframes 
data_day_year <- data_day_allyears
data_n_year   <- data_n_allyears
#rename dataframes 
library(xtable)
library(dplyr)
library(rmeta)
library(rdrobust)

year <- c(2009, 2013, 2014, 2015)

data <- left_join(data_day_year, data_n_year,
                  by = c("days", "year"))

data_close <- 
  data %>%
  mutate(voted     = round(par_vote*n),
         abstained = round((1-par_vote)*n),
         treated   = days > 0) %>% 
  select(year, days, treated, voted, abstained)

data_voted <- 
  as.data.frame(lapply(data_close, 
                       function(x,p) rep(x,p), 
                       data_close[["voted"]])) %>%
  mutate(voted = 1)

data_abstained <- 
  as.data.frame(lapply(data_close, 
                       function(x,p) rep(x,p), 
                       data_close[["abstained"]])) %>%
  mutate(voted = 0)

data <- 
  rbind(data_voted, data_abstained)

# array with coefficients for years and models 

main_table_out <- matrix(NA, nrow = 11, ncol = 5)

coef_matrix <- matrix(NA, nrow = 4, ncol = 5)
for (i in 1:4){
  year_vec <- year[i]
  data_year <- 
    data %>%
    ungroup() %>%
    mutate(days = days) %>% 
    filter(year == year_vec )
  
  rd <- rdrobust(data_year$voted, data_year$days)
  coef_matrix[c(1,3), i] <- rd$coef[c(1,3)]
  coef_matrix[c(2,4), i] <- rd$se[c(1,3)]
  
  main_table_out[7, i] <- as.numeric(data_year %>% 
                                       filter(treated == 0 & abs(days) < rd$h_l) %>% 
                                       summarise(mean(voted)))
  main_table_out[8  ,i] <- rd$N_h_l
  main_table_out[9  ,i] <- rd$N_h_r
  main_table_out[10    ,i] <- round(rd$h_l)
  main_table_out[11    ,i] <- round(rd$b_l)
  print(i)
}

# Pool coefficients using fixed effects meta analysis
coef_matrix[1:2, 5] <- c(meta.summaries(coef_matrix[ 1, 1:4],
                                        coef_matrix[ 2, 1:4])$summary,
                         meta.summaries(coef_matrix[ 1, 1:4],
                                        coef_matrix[ 2, 1:4])$se.summary)
coef_matrix[3:4, 5] <- c(meta.summaries(coef_matrix[ 3, 1:4],
                                        coef_matrix[ 4, 1:4])$summary,
                         meta.summaries(coef_matrix[ 3, 1:4],
                                        coef_matrix[ 4, 1:4])$se.summary)

main_table_out[7,5] <- weighted.mean(x = main_table_out[7, 1:4],
                                     w = main_table_out[8, 1:4])
main_table_out[7,] <- round(100*main_table_out[7,], 2)
main_table_out[c(8,9),5] <- rowSums(main_table_out[c(8,9), -5])
for (j in 1:5){
  main_table_out[1,j] <- 100*round(coef_matrix[1,j], 4)
  main_table_out[2,j] <- 100*round(coef_matrix[2,j], 4)
  main_table_out[4,j] <- 100*round(coef_matrix[3,j], 4)
  main_table_out[5,j] <- 100*round(coef_matrix[4,j], 4)
  main_table_out[3     , j] <- paste("[",
                                     100*round(coef_matrix[1,j] + coef_matrix[2,j] * qnorm(0.025), 4),
                                     ", ",
                                     100*round(coef_matrix[1,j] + coef_matrix[2,j] * qnorm(0.975), 4),
                                     "]",
                                     sep = "")
  main_table_out[6     , j] <- paste("[",
                                     100*round(coef_matrix[3,j] + coef_matrix[4,j] * qnorm(0.025), 4),
                                     ", ",
                                     100*round(coef_matrix[3,j] + coef_matrix[4,j] * qnorm(0.975), 4),
                                     "]",
                                     sep = "")  
  
}

rownames(main_table_out) <- c("Estimate", "St.err.", "CI",
                              "Robust estimate", "Robust st.err.", "Robust CI",
                              "Mean turnout among parents of ineligible", "N_ineligible",
                              "N_eligible", "Bandwidth", "Bias-correction bandwidth")
colnames(main_table_out) <- c("2009", "2013", "2014", "2015", "Pooled Effect")

main_table_out <- rbind(main_table_out_all, main_table_out)

xtable(main_table_out)

